function [a,b,R2] = generic_nominal_reg(inflation_process, lhs0, t, lhst, rhs0, prob_NxYz, Phi_NxYz_t, grid_Y)
% regression is of the following form:
% lhs(0) + lhs(t) = a + b*rhs(0) + err
% where the terms F(t)=lhs(0), lhs(t), rhs(0) are all of the form
% F(t) = fun_Nxz(N(t),x(t),z(t)) + const + slope_Y*Y(t) + slope_X*X(t) + slope_eps*eps(t)

    if t<=0; error('Expecting t>0'); end

    [NN, Nx, NY, Nz] = size(prob_NxYz);

    mean_X = inflation_process.mu_pi/(1 - inflation_process.rho_pi);
    var_X = (1 + 2*inflation_process.rho_pi*inflation_process.theta_pi ...
        + inflation_process.theta_pi^2)*(inflation_process.sigma_pi^2)...
        /(1 - inflation_process.rho_pi^2);
    var_eps = inflation_process.sigma_pi^2;
    cov_Xe = var_eps;
    
    cov_X0_Xt = (inflation_process.rho_pi^t)*var_X ...
        + (inflation_process.rho_pi^(t-1))*inflation_process.theta_pi*cov_Xe;
    cov_eps0_Xt = (inflation_process.rho_pi^t)*cov_Xe ...
        + (inflation_process.rho_pi^(t-1))*inflation_process.theta_pi*var_eps;
    
    Y_NxYz = permute(repmat(grid_Y(:),[1 NN Nx Nz]),[2 3 1 4]);
    
    % 1. rhs
    rhs_NxYz = permute(repmat(rhs0.fun_Nxz,[1 1 1 NY]),[1 2 4 3]) ...
        + rhs0.slope_Y*Y_NxYz;
    mean_rhs_NxYz = sum(prob_NxYz.*rhs_NxYz,'all');
    var_rhs_NxYz = sum(prob_NxYz.*(rhs_NxYz.^2),'all') - mean_rhs_NxYz^2;
    
    mean_rhs_Xe = rhs0.const + rhs0.slope_X*mean_X;
    var_rhs_Xe = (rhs0.slope_X^2)*var_X + (rhs0.slope_eps^2)*var_eps ...
        + 2*rhs0.slope_X*rhs0.slope_eps*cov_Xe;
    
    mean_rhs = mean_rhs_NxYz + mean_rhs_Xe;
    var_rhs = var_rhs_NxYz + var_rhs_Xe;
    
    % 2. lhs0
    lhs0_NxYz = permute(repmat(lhs0.fun_Nxz,[1 1 1 NY]),[1 2 4 3]) ...
        + lhs0.slope_Y*Y_NxYz;
    mean_lhs0_NxYz = sum(prob_NxYz.*lhs0_NxYz,'all');
    var_lhs0_NxYz = sum(prob_NxYz.*(lhs0_NxYz.^2),'all') - mean_lhs0_NxYz^2;
    
    mean_lhs0_Xe = lhs0.const+ lhs0.slope_X*mean_X;
    var_lhs0_Xe = (lhs0.slope_X^2)*var_X + (lhs0.slope_eps^2)*var_eps ...
        + 2*lhs0.slope_X*lhs0.slope_eps*cov_Xe;
    
    % 3. lhst
    lhst_NxYz = permute(repmat(lhst.fun_Nxz,[1 1 1 NY]),[1 2 4 3]) ...
        + lhst.slope_Y*Y_NxYz;
    mean_lhst_NxYz = sum(prob_NxYz.*lhst_NxYz,'all');
    var_lhst_NxYz = sum(prob_NxYz.*(lhst_NxYz.^2),'all') - mean_lhst_NxYz^2;
    
    mean_lhst_Xe = lhst.const+ lhst.slope_X*mean_X;
    var_lhst_Xe = (lhst.slope_X^2)*var_X + (lhst.slope_eps^2)*var_eps ...
        + 2*lhst.slope_X*lhst.slope_eps*cov_Xe;
    
    % 4. lhs
    cov_lhs0_lhst_NxYz = sum(prob_NxYz(:).*lhs0_NxYz(:).*(Phi_NxYz_t*lhst_NxYz(:))) ...
        - mean_lhs0_NxYz*mean_lhst_NxYz;
    cov_lhs0_lhst_Xe = lhs0.slope_X*lhst.slope_X*cov_X0_Xt ...
        + lhs0.slope_eps*lhst.slope_X*cov_eps0_Xt;
    var_lhs = var_lhs0_NxYz + var_lhs0_Xe + var_lhst_NxYz + var_lhst_Xe ...
        + 2*cov_lhs0_lhst_NxYz + 2*cov_lhs0_lhst_Xe;
    
    % 5. cov(lhs,rhs)
    cov_lhs0_rhs_NxYz = sum(prob_NxYz(:).*lhs0_NxYz(:).*rhs_NxYz(:)) ...
        - mean_lhs0_NxYz*mean_rhs_NxYz;
    cov_lhst_rhs_NxYz = sum(prob_NxYz(:).*rhs_NxYz(:).*(Phi_NxYz_t*lhst_NxYz(:))) ...
        - mean_lhst_NxYz*mean_rhs_NxYz;
    
    cov_lhs0_rhs_Xe = lhs0.slope_X*rhs0.slope_X*var_X + (lhs0.slope_X*rhs0.slope_eps ...
        + lhs0.slope_eps*rhs0.slope_X)*cov_Xe + lhs0.slope_eps*rhs0.slope_eps*var_eps;
    cov_lhst_rhs_Xe = lhst.slope_X*rhs0.slope_X*cov_X0_Xt + lhst.slope_X*rhs0.slope_eps*cov_eps0_Xt;
    
    cov_lhs_rhs = cov_lhs0_rhs_NxYz + cov_lhst_rhs_NxYz + cov_lhs0_rhs_Xe ...
        + cov_lhst_rhs_Xe;
    
    % 6. reg coeffs
    b = cov_lhs_rhs/var_rhs;
    a = mean_lhs0_NxYz + mean_lhst_NxYz + mean_lhs0_Xe + mean_lhst_Xe - b*mean_rhs;
    R2 = (b^2*var_rhs)/var_lhs;
    

end